---
title: "CARRS RDD analysis using a standardized cutoff - pair wise results"
knit: (function(input, ...) {rmarkdown::render(input, output_file="")})
output: pdf_document
geometry: "left=2cm,right=2cm,top=1cm,bottom=1cm"
date: "`r format(Sys.time(), '%d %B, %Y')`"
---


```{r , setup, include=FALSE, eval = TRUE}
knitr::opts_chunk$set(
   comment = '', warning=FALSE,  message=FALSE, results= FALSE , fig.width = 9,fig.height = 7, echo=FALSE, root.dir =""
)

 knitr::opts_knit$set(
   root.dir = ""    
 )
```

```{r packages}
rm(list=ls())
  library(tidyverse)
  library(gridExtra)
  library(dplyr)
  library(rdrobust)
  library(rddensity)
  library(viridisLite)
  library(foreign)
  library(ggplot2)
  library(data.table)
  library(doBy)
  library(socviz)
  library(openxlsx) 
  library(modelsummary)
  library(labelled)
  library(gtsummary)
  library(gt)
  library(kableExtra)
```

```{r loaddata}
data_treatment<-data.frame(read.csv("carrs_full.csv"))
data_treatment<-as.data.table(data_treatment)

data_diagnosis<-data.frame(read.csv("carrs_selected.csv"))
data_diagnosis<-as.data.table(data_diagnosis)
```


##calculating the bandwidth for each sub group
```{r}
# Males and females separately
rv.update = function(a,b,c) {
  out=c[a==1 & female==b,runningvar:=ifelse((((sys-140)/(sd(sys))) == ((dias-90)/(sd(dias)))),
                                             (sys-140)/(sd(sys)),
                                            ifelse((((sys-140)/(sd(sys))) > ((dias-90)/(sd(dias)))),((sys-140)/(sd(sys))),((dias-90)/(sd(dias)))))]
a_b<<-out
}

rv.update('w01',1,data_diagnosis)
w01_female_diag<-a_b[w01==1 & female==1,]
rm(a_b)

rv.update('w01',0,data_diagnosis)
w01_male_diag<-a_b[w01==1 & female==0,]
rm(a_b)

rv.update('w23',1,data_diagnosis)
w23_female_diag<-a_b[w23==1 & female==1,]
rm(a_b)

rv.update('w23',0,data_diagnosis)
w23_male_diag<-a_b[w23==1 & female==0,]
rm(a_b)

rv.update('w45',1,data_diagnosis)
w45_female_diag<-a_b[w45==1 & female==1,]
rm(a_b)

rv.update('w45',0,data_diagnosis)
w45_male_diag<-a_b[w45==1 & female==0,]
rm(a_b)

rv.update('w01',1,data_treatment)
w01_female_treat<-a_b[w01==1 & female==1,]
rm(a_b)

rv.update('w01',0,data_treatment)
w01_male_treat<-a_b[w01==1 & female==0,]
rm(a_b)

rv.update('w23',1,data_treatment)
w23_female_treat<-a_b[w23==1 & female==1,]
rm(a_b)

rv.update('w23',0,data_treatment)
w23_male_treat<-a_b[w23==1 & female==0,]
rm(a_b)

rv.update('w45',1,data_treatment)
w45_female_treat<-a_b[w45==1 & female==1,]
rm(a_b)

rv.update('w45',0,data_treatment)
w45_male_treat<-a_b[w45==1 & female==0,]
rm(a_b)


# Total population
rv.updateMT = function(a,c) {
  out=c[a==1,runningvar:=ifelse((((sys-140)/(sd(sys))) == ((dias-90)/(sd(dias)))),
                                             (sys-140)/(sd(sys)),
                                            ifelse((((sys-140)/(sd(sys))) > ((dias-90)/(sd(dias)))),((sys-140)/(sd(sys))),((dias-90)/(sd(dias)))))]
  a_b<<-out
}

# diagnosis
rv.updateMT('w01',data_diagnosis)
w01_total_diag<-a_b[w01==1 ,]
rm(a_b)

rv.updateMT('w23',data_diagnosis)
w23_total_diag<-a_b[w23==1 ,]

rm(a_b)

rv.updateMT('w45',data_diagnosis)
w45_total_diag<-a_b[w45==1 ,]
rm(a_b)

# treatment
rv.updateMT('w01',data_treatment)
w01_total_treat<-a_b[w01==1,]
rm(a_b)

rv.updateMT('w23',data_treatment)
w23_total_treat<-a_b[w23==1,]
rm(a_b)


rv.updateMT('w45',data_treatment)
w45_total_treat<-a_b[w45==1,]
rm(a_b)
```
##RDD Outputs function

```{r fuc}
  carrs_rdd = function(a,b,d) {
    d$Y<-a
    d$X<-b

    out =  rdrobust(y = d$Y, d$X, kernel = "triangular", p = 1, 
                    c = 0, rho = 1 , all = T)
    summary(out)
    
    ret <- data.frame(
sex = c(formatC(round(100*out$coef[1,1],2),2,format="f"), formatC(round(out$pv[3,1],2),2,format="f"), 
            paste0( "(", formatC(round(100*out$ci[3,1],2),2,format="f"), " - ", formatC(round(100*out$ci[3,2],2),2,format="f"), ")" ), 
            formatC(round(out$bws[1,1],2),2,format="f"), formatC(out$N_h[1] + out$N_h[2]))
  )
    rownames(ret) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
    
    output<<-ret %>% add_rownames  
  }
```

##Calling the specific subroups and saving the output

```{r}
rdd = function(a) {
  carrs_rdd(a$hypdiag,a$runningvar,a)
}

rdd(w01_male_diag)
w01_male_diag <-output %>% mutate(Male=sex)%>%select(-sex)

rdd(w23_male_diag)
w23_male_diag <-output%>%mutate(Male=sex)%>%select(-sex)

rdd(w45_male_diag)
w45_male_diag <-output%>%mutate(Male=sex)%>%select(-sex)


rdd(w01_female_diag)
w01_female_diag <-output%>%mutate(Female=sex)%>%select(-sex)

rdd(w23_female_diag)
w23_female_diag <-output%>%mutate(Female=sex)%>%select(-sex)

rdd(w45_female_diag)
w45_female_diag <-output%>%mutate(Female=sex)%>%select(-sex)


rdd(w01_total_diag)
w01_total_diag <-output%>%mutate(Total=sex)%>%select(-sex)

rdd(w23_total_diag)
w23_total_diag <-output%>%mutate(Total=sex)%>%select(-sex)

rdd(w45_total_diag)
w45_total_diag <-output%>%mutate(Total=sex)%>%select(-sex)
```

```{r treatment}
rdd = function(a) {
  carrs_rdd(a$hypdrug,a$runningvar,a)
}

rdd(w01_male_treat)
w01_male_treat <-output%>%mutate(Male=sex)%>%select(-sex)

rdd(w23_male_treat)
w23_male_treat <-output%>%mutate(Male=sex)%>%select(-sex)

rdd(w45_male_treat)
w45_male_treat <-output%>%mutate(Male=sex)%>%select(-sex)


rdd(w01_female_treat)
w01_female_treat <-output%>%mutate(Female=sex)%>%select(-sex)

rdd(w23_female_treat)
w23_female_treat <-output%>%mutate(Female=sex)%>%select(-sex)

rdd(w45_female_treat)
w45_female_treat <-output%>%mutate(Female=sex)%>%select(-sex)


rdd(w01_total_treat)
w01_total_treat <-output%>%mutate(Total=sex)%>%select(-sex)

rdd(w23_total_treat)
w23_total_treat <-output%>%mutate(Total=sex)%>%select(-sex)

rdd(w45_total_treat)
w45_total_treat <-output%>%mutate(Total=sex)%>%select(-sex)
```

```{r}
table_results = function(x,y,z) {
  x<-x
  y<-y
  z<-z
  c<-merge(x,y,by="rowname")
  c<-merge(c,z,by="rowname")
  c<-c%>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))
out<<-c
}
```



```{r}
table_results(w01_total_diag,w23_total_diag,w45_total_diag) 
pairs_table_tdiag<-out
rm(c)
colnames(pairs_table_tdiag)[2] <- "Pair 1" 
colnames(pairs_table_tdiag)[3] <- "Pair 2"   
colnames(pairs_table_tdiag)[4] <- "Pair 3"   
colnames(pairs_table_tdiag)[1] <- " " 


table_results(w01_female_diag,w23_female_diag,w45_female_diag)
pairs_table_fdiag<-out
rm(c)
colnames(pairs_table_fdiag)[2] <- "Pair 1" 
colnames(pairs_table_fdiag)[3] <- "Pair 2"   
colnames(pairs_table_fdiag)[4] <- "Pair 3"   
colnames(pairs_table_fdiag)[1] <- " " 


table_results(w01_male_diag,w23_male_diag,w45_male_diag)
pairs_table_mdiag<-out
rm(c)
colnames(pairs_table_mdiag)[1] <- " " 
colnames(pairs_table_mdiag)[2] <- "Pair 1" 
colnames(pairs_table_mdiag)[3] <- "Pair 2"   
colnames(pairs_table_mdiag)[4] <- "Pair 3"   

comb_diag_table<-rbind(pairs_table_tdiag,pairs_table_fdiag,pairs_table_mdiag)


table_results(w01_total_treat,w23_total_treat,w45_total_treat) 
pairs_table_ttreat<-out
rm(c)
colnames(pairs_table_ttreat)[2] <- "Pair 1" 
colnames(pairs_table_ttreat)[3] <- "Pair 2"   
colnames(pairs_table_ttreat)[4] <- "Pair 3"   
colnames(pairs_table_ttreat)[1] <- " " 


table_results(w01_female_treat,w23_female_treat,w45_female_treat)
pairs_table_ftreat<-out
rm(c)
colnames(pairs_table_ftreat)[2] <- "Pair 1" 
colnames(pairs_table_ftreat)[3] <- "Pair 2"   
colnames(pairs_table_ftreat)[4] <- "Pair 3"   
colnames(pairs_table_ftreat)[1] <- " " 


table_results(w01_male_treat,w23_male_treat,w45_male_treat)
pairs_table_mtreat<-out
rm(c)
colnames(pairs_table_mtreat)[1] <- " " 
colnames(pairs_table_mtreat)[2] <- "Pair 1" 
colnames(pairs_table_mtreat)[3] <- "Pair 2"   
colnames(pairs_table_mtreat)[4] <- "Pair 3"   

comb_treat_table<-rbind(pairs_table_ttreat,pairs_table_ftreat,pairs_table_mtreat)

comb_pair_table<-cbind(comb_diag_table,comb_treat_table)
comb_pair_table<-comb_pair_table[ , -c(5) ]

colnames(comb_pair_table) <-c('','Pair 1','Pair 2', 'Pair 3',
                                           'Pair 1','Pair 2', 'Pair 3')
```

#combining diagnosis and treatment
```{r }
  filename <- paste0("pairwise_combined.tex")
  caption=paste0("Effect of home-based screening on hypertension diagnosis and treatment by survey pair")

x = kable(head(mtcars), "html")
combined<-kbl(comb_pair_table,caption=caption, booktabs = T, linesep = '', escape = FALSE,align = "lccccccccc", format = "latex",digits = 2) %>%
    add_header_above(c(" "=1, "Diagnosis" = 3, "Treatment" = 3)) %>%
    kable_styling(latex_options = "HOLD_position", position = "center") %>%
    pack_rows("Panel A: Total", 1,5) %>%
    pack_rows("Estimation results", 1,3) %>%
    pack_rows("Bandwidth details", 4, 5) %>%
    pack_rows(" ", 6,10) %>%
    pack_rows("Panel B: Female", 6,10) %>%
    pack_rows("Estimation results", 6, 8) %>%
    pack_rows("Bandwidth details", 9, 10) %>%
    pack_rows(" ", 11,15) %>%
    pack_rows("Panel C: Male", 11,15) %>%
    pack_rows("Estimation results", 11, 13) %>%
    pack_rows("Bandwidth details", 14, 15) %>%
    footnote(general = c("The model estimation included a local linear trend, triangular Kernel weights, and a data-driven bandwidth selection (mean squared error optimal bandwidth). The table displays the local average treatment effect in percentage points for individuals with a standardized blood pressure of zero and CIs resulting from robust bias-corrected standard errors clustered at the individual-level.","Abbreviation: CI = confidence interval"),general_title="",threeparttable = T, fixed_small_size = T)%>%
    save_kable(filename, format = "latex")
```

